#include <iostream>
#include <fstream>

#include "TH1.h"
#include "TF1.h"
#include "TFile.h"
#include "TTree.h"
#include "Math/WrappedTF1.h"
#include "Math/GSLIntegrator.h"
#include "TMath.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TImage.h"
#include "TVector3.h"
#include "TH3.h"

using namespace std;

int main(int argc, char** argv)
{
	gROOT->Reset();
	gROOT->SetStyle("Plain"); // Start with clean sheet for plotting graphs.

	//
	// Increment distance in cm of plotting vector for fermat surface
	// reconstruction.
	//
	const Double_t plotVectorIncrement = 400.;
	//
	// Effective index of refraction in liquid
	//scintillator. This value should be ~ 1.44.
	//
	const Double_t index = 1.44;
	const Double_t speedOfLight = 2.99792458e8*1.e+2*1.e-9; // Speed of light in vacuum in cm per ns.
	//
	// Radius of Detector in cm. Actually the real radius of the detector is 900
	// cm but 850 cm is the radius where the largest slice-section of the PMT is
	// placed WRT to center of the volume.
	//
	const Double_t radiusOfDetector = 850.; // Radius of Detector in cm.
	const Double_t speedOfLightInLS = speedOfLight/index;
	//	const Double_t transitTimeSpread = 50.; // Average TTS of PMT in ns.
	//	const Double_t transitTimeSpread = 5.; // Average TTS of PMT in ns.
	const Int_t maxSize = 10000;

	TVector3 trackBeginningVector;
	Double_t trackBeginningVectorFOM;
	TVector3 trackEndVector;
	TTree* plotTree = new TTree("plotTree", "plot data");
	Double_t plotVectorX, plotVectorY, plotVectorZ;
	Double_t FOM, maxFOM, minFOM;
	plotTree->Branch("X component of plot vector", &plotVectorX, "plotVectorX/D");
	plotTree->Branch("Y component of plot vector", &plotVectorY, "plotVectorY/D");
	plotTree->Branch("Z component of plot vector", &plotVectorZ, "plotVectorZ/D");
	plotTree->Branch("Emission likelyhood Figure of Merit", &FOM, "FOM/D");
	Double_t reconThetaFromRealAngle;
	TTree* reconstructedTrackVectorDivergenceAngleTree = new TTree("reconThetaFromRealAngle", "divergance angles");
	reconstructedTrackVectorDivergenceAngleTree->Branch("reconThetaFromRealAngle", &reconThetaFromRealAngle, "reconThetaFromRealAngle/D");
	Int_t PMTCableNo[maxSize];
	Double_t PMTHitTime[maxSize];
	//	Double_t PMTCharge[maxSize];
	TVector3 PMTPositionVector[maxSize];
	TVector3 plotVector(0., 0., 0.);
	Int_t cableNo, surfaceNumber, PMTiNumber, PMTjNumber;
	Double_t x, y, z, t;
	for(Int_t i = 0; i <= maxSize; i++){ // Initialize arrays.
		PMTCableNo[i] = -1;
		// Set some negative number as flag.
		PMTHitTime[i] = TMath::SignalingNaN();
		PMTPositionVector[i].SetXYZ(TMath::SignalingNaN(), TMath::SignalingNaN(), TMath::SignalingNaN());
	}

	cout << endl;
	cout << endl;
	cout << " **************************************" << endl;
	cout << " *                                    *" << endl;
	cout << " * Event Image Reconstruction Program *" << endl;
	cout << " *                                    *" << endl;
	cout << " **************************************" << endl;
	cout << endl;
	cout << endl;
	const char* inputDataFile = argv[1];
	ifstream serialFIn(inputDataFile); // input file holding pmt position and timing info
	char charBuffer[256];
	string bufferString;
	while(serialFIn.getline(charBuffer, 256)){
		bufferString = charBuffer;
		if(bufferString[bufferString.find_first_not_of("	 ")] == '#') continue;
		istringstream streamIn(charBuffer);
		if(streamIn >> cableNo >> cableNo >> surfaceNumber >> PMTiNumber >> PMTjNumber >> x >> y >> z >> t){
			//PMTCableNo[cableNo] = cableNo;
			PMTPositionVector[cableNo].SetXYZ(x, y, z);
			//PMTHitTime[cableNo] = t;
		}
		else{
			cout << "Error: file format is not correct." << endl;
			cout << "Exiting." << endl;
			return 0;
		}
	}
	serialFIn.close();

	//Some variableT holders
	Int_t run;
	Int_t eventNumber;
	Double_t timeStamp;
	Float_t mpArr[3];//3 array of xyz of track midpoint
	Float_t dirnArr[3];//3 array of xyz of track direction
	Int_t nPulses;//there are 4 branches which are arrays of size N=nPulses
	Int_t cableArr[maxSize];// pointer for cable number array
	Int_t gainArr[maxSize];//pointer for gain array
	Int_t AorBArr[maxSize];//pointer for AorB array;
	Double_t hitTimeArr[maxSize];//pointer for array of PMT hit times
	Double_t hitQArr[maxSize];//pointer for array of PMT hit times

	TFile* dataFIn = new TFile("testDataGrabber2.root","READ");
	if(!Bool_t(dataFIn) || !dataFIn->IsOpen()){
		cout << "Oops, there is a problem with the file, try again ! " << endl;
		return 0;
	}
	TTree* dataTree  = (TTree*)dataFIn->Get("KLMuonsAndPMTHitTimes");
	if(!dataTree){
		cout << "Oops, there is a problem with the tree, please check it! " << endl;
		return 0;
	}
	dataTree->SetBranchAddress("KLRun",&run);
	dataTree->SetBranchAddress("KLEventNumber",&eventNumber);
	dataTree->SetBranchAddress("KLTimestamp",&timeStamp);
	dataTree->SetBranchAddress("KLTrackMidpoint",mpArr);
	dataTree->SetBranchAddress("KLTrackDirection",dirnArr);
	dataTree->SetBranchAddress("NPulses",&nPulses);
	dataTree->SetBranchAddress("CableNum", cableArr);
	dataTree->SetBranchAddress("ATWDAorB", AorBArr);
	dataTree->SetBranchAddress("channelGain", gainArr);
	dataTree->SetBranchAddress("PMTHitTime", hitTimeArr);
	dataTree->SetBranchAddress("PMTHitQ", hitQArr);

	//		cout<<"Here is the data for entry "<<i<<endl;
	//		cout<<"KamLAND Run = "<<run<<endl;
	//		cout<<"EventNumber = "<<eventNumber<<endl;
	//		cout<<"Timestamp = "<<timeStamp<<endl;
	//		cout<<"nPulses = "<<nPulses<<endl;
	//		//Track midpoint
	//		cout<<"The KL track midpoint is "<<endl;
	//		cout<<"X: "<<mpArr[0]<<endl;
	//		cout<<"Y: "<<mpArr[1]<<endl;
	//		cout<<"Z: "<<mpArr[2]<<endl;
	//		//Track direction     
	//		cout<<"The KL track direction is "<<endl;
	//		cout<<"X: "<<dirnArr[0]<<endl;
	//		cout<<"Y: "<<dirnArr[1]<<endl;
	//		cout<<"Z: "<<dirnArr[2]<<endl;
	//		cout<<"Entry\tCable\tGain\tAorB\tHitTime\thitCharge(npe)"<<endl;
	//		for(Int_t j=0;j<=nPulses;j++){
	//			if(j>=nPulses)continue;//make sure array index is less than size of array, nPulses
	//			cout<<j<<"\t"<<cableArr[j]<<"\t"<<gainArr[j]<<"\t"<<AorBArr[j]<<"\t"<<hitTimeArr[j]<<"\t"<<hitQArr[j]<<endl;
	//		}

	// One sigma value. Here we probably want the time for
	// one photon to travel between ajacent PMTs.
	TVector3 differenceVector;
	Double_t deltaT, deltaTMin, deltaTMax;
	const Double_t sigma = 1.;
	Double_t tIntegrationIncrement = sigma;
	Double_t function;
	Int_t numerator = radiusOfDetector/plotVectorIncrement;
	Double_t plotVectorStartingPosition = -(numerator)*plotVectorIncrement;

	for(Int_t i = 0; i < dataTree->GetEntries(); i++){
		maxFOM = TMath::SignalingNaN();
		minFOM = TMath::SignalingNaN();
		trackBeginningVectorFOM = TMath::SignalingNaN();
		trackBeginningVector.SetXYZ(TMath::SignalingNaN(), TMath::SignalingNaN(), TMath::SignalingNaN());
		dataTree->GetEntry(i);
		cout
			<< endl << endl
			<< "Starting analysis for event "
			<< i
			<< ".\n";
		for(Int_t j = 0; j < maxSize; j++){
			PMTHitTime[j] = TMath::SignalingNaN();
		}
		for(Int_t j = 0; j < nPulses; j++){
			PMTHitTime[cableArr[j]] = hitTimeArr[j];
		}
		for(plotVector.SetX(plotVectorStartingPosition); plotVector.X() <= radiusOfDetector; plotVector.SetX(plotVector.X() + plotVectorIncrement)){
			//	for(plotVector.SetX(0.); plotVector.X() == 0.; plotVector.SetX(plotVector.X() + plotVectorIncrement)){
			for(plotVector.SetY(plotVectorStartingPosition); plotVector.Y() <= radiusOfDetector; plotVector.SetY(plotVector.Y() + plotVectorIncrement)){
				//		for(plotVector.SetY(0.); plotVector.Y() == 0.; plotVector.SetY(plotVector.Y() + plotVectorIncrement)){
				for(plotVector.SetZ(plotVectorStartingPosition); plotVector.Z() <= radiusOfDetector; plotVector.SetZ(plotVector.Z() + plotVectorIncrement)){
					if(plotVector.Mag() > radiusOfDetector) continue;
					function = 0.;
					FOM = 0.;
					deltaTMin = TMath::SignalingNaN();
					deltaTMax = TMath::SignalingNaN();
					for(Int_t j = 0; j < maxSize; j++){
						if(TMath::IsNaN(PMTHitTime[j])) continue;
						differenceVector = PMTPositionVector[j] - plotVector;
						deltaT = PMTHitTime[j] -	(differenceVector.Mag())/speedOfLightInLS;
						if(TMath::IsNaN(deltaTMax)) deltaTMax = deltaT;
						else if(deltaTMax < deltaT) deltaTMax = deltaT;
						if(TMath::IsNaN(deltaTMin)) deltaTMin = deltaT;
						else if(deltaTMin > deltaT) deltaTMin = deltaT;
					}
//					cout
//						<< "deltaTMin = " << deltaTMin << endl;
//					cout
//						<< "deltaTMax = " << deltaTMax << endl;
					for(Double_t variableT = deltaTMin - 5.*sigma; variableT <= deltaTMax + 5.*sigma; variableT += tIntegrationIncrement){
						function = 0.;
						for(Int_t j = 0; j <= maxSize; j++){
							if(TMath::IsNaN(PMTHitTime[j])) continue;
							differenceVector = PMTPositionVector[j] - plotVector;
							deltaT = PMTHitTime[j] -	(differenceVector.Mag())/speedOfLightInLS;
							function +=
								(variableT - deltaT)
								*exp(-pow(variableT	- deltaT, 2)
										/(2.*pow(sigma, 2)));
						}
						FOM += tIntegrationIncrement * pow(function, 2);
					}
					if(TMath::IsNaN(trackBeginningVectorFOM)){
						trackBeginningVectorFOM = FOM;
						trackBeginningVector = plotVector;
					}
					else if(trackBeginningVectorFOM < FOM){
						trackBeginningVectorFOM = FOM;
						trackBeginningVector = plotVector;
					}
					if(TMath::IsNaN(minFOM)) minFOM = FOM;
					else if(minFOM > FOM) minFOM = FOM;
					if(TMath::IsNaN(maxFOM)) maxFOM = FOM;
					else if(maxFOM < FOM) maxFOM = FOM;

					plotVectorX = plotVector.X();
					plotVectorY = plotVector.Y();
					plotVectorZ = plotVector.Z();
					plotTree->Fill();
					cout
						<< "XYZ coord: "
						<< plotVector.X()
						<< ", "
						<< plotVector.Y()
						<< ", "
						<< plotVector.Z()
						<< " => Photon emission FOM: "
						<< FOM
						<< endl;
				}
			}
		}

		cout
			<< "Most likely event beginning vertex: "
			<< trackBeginningVector.X()
			<< ", "
			<< trackBeginningVector.Y()
			<< ", "
			<< trackBeginningVector.Z()
			<< endl;
		
		Double_t FOMThreshold = minFOM + .2*(maxFOM - minFOM)/*1500.*/;
		cout
			<< "Minimum FOM: "
			<< minFOM
			<< ", "
			<< "Maximum FOM: "
			<< maxFOM
			<< ", "
			<< "FOM threshold: "
			<< FOMThreshold
			<< endl;
		cout
			<< "KL track midpoint = ("
			<< mpArr[0] << ", "
			<< mpArr[1] << ", "
			<< mpArr[2] << ")\n";
		cout
			<< "KL track direction = ("
			<< dirnArr[0] << ", "
			<< dirnArr[1] << ", "
			<< dirnArr[2] << ")\n";

		//
		// Reconstructed track fitting method 1.
		//
		//TVector3 testTrackEndVector, testTrackVector, variableVector,
		//		 deviationVector;
		//Double_t deviation = TMath::SignalingNaN(), leastDeviation = TMath::SignalingNaN();
		//for(Int_t j = 0; j < plotTree->GetEntries(); j++){
		//	plotTree->GetEntry(j);
		//	if(FOM < FOMThreshold) continue;
		//	testTrackEndVector.SetXYZ(plotVectorX, plotVectorY, plotVectorZ);
		//	testTrackVector = testTrackEndVector - trackBeginningVector;
		//	deviation = 0.;
		//	for(Int_t k = 0; k < plotTree->GetEntries(); k++){
		//		plotTree->GetEntry(k);
		//		if(FOM < FOMThreshold) continue;
		//		variableVector.SetXYZ(plotVectorX, plotVectorY, plotVectorZ);
		//		deviationVector = variableVector - trackBeginningVector;
		//		Double_t projectionOfVariableVector =
		//			deviationVector.Dot(testTrackVector)/testTrackVector.Mag();
		//		if(projectionOfVariableVector < 0.){
		//			deviation +=
		//				FOM*deviationVector.Mag();
		//		}
		//		else if(projectionOfVariableVector >= 0. &&
		//				projectionOfVariableVector <=
		//				testTrackVector.Mag()){
		//			deviationVector = (1./testTrackVector.Mag())*(testTrackVector.Cross(deviationVector));
		//			deviation +=
		//				FOM*deviationVector.Mag();
		//		}
		//		else{
		//			deviationVector = deviationVector - testTrackVector;
		//			deviation +=
		//				FOM*deviationVector.Mag();
		//		}
		//	}
		//	if(TMath::IsNaN(leastDeviation)){
		//		leastDeviation = deviation;
		//		trackEndVector = testTrackEndVector;
		//	}
		//	else if(leastDeviation > deviation){
		//		trackEndVector = testTrackEndVector;
		//		leastDeviation = deviation;
		//	}
		//}

		//Double_t fittingVectorIncrement = plotVectorIncrement;
		Double_t deviationMin = TMath::SignalingNaN();
		//
		// Reconstructed track fitting method 2.
		//
		TVector3 testReconTrackBeginningVector,
				 testReconTrackEndVector,
				 reconstructedTrackBeginningVector,
				 reconstructedTrackEndVector;
		//for(Int_t j = 0; j < plotTree->GetEntries(); j++){
		//	plotTree->GetEntry(j);
		//	if(FOM < FOMThreshold) continue;
		//	testReconTrackBeginningVector.SetXYZ(plotVectorX, plotVectorY, plotVectorZ);
		//	for(Int_t k = 0; k < plotTree->GetEntries(); k++){
		//		plotTree->GetEntry(k);
		//		if(FOM < FOMThreshold) continue;
		//		testReconTrackEndVector.SetXYZ(plotVectorX, plotVectorY, plotVectorZ);
		//		if(testReconTrackEndVector == testReconTrackBeginningVector) continue;
		//numerator = radiusOfDetector/fittingVectorIncrement;
		//Double_t testReconTrackBeginningVectorStartingPosition = -(numerator)*fittingVectorIncrement,
		//		 testReconTrackEndVectorStartingPosition = -(numerator)*fittingVectorIncrement;
		//for(testReconTrackBeginningVector.SetX(testReconTrackBeginningVectorStartingPosition); testReconTrackBeginningVector.X() <= radiusOfDetector; testReconTrackBeginningVector.SetX(testReconTrackBeginningVector.X() + fittingVectorIncrement)){
		//	for(testReconTrackBeginningVector.SetY(testReconTrackBeginningVectorStartingPosition); testReconTrackBeginningVector.Y() <= radiusOfDetector; testReconTrackBeginningVector.SetY(testReconTrackBeginningVector.Y() + fittingVectorIncrement)){
		//		for(testReconTrackBeginningVector.SetZ(testReconTrackBeginningVectorStartingPosition); testReconTrackBeginningVector.Z() <= radiusOfDetector; testReconTrackBeginningVector.SetZ(testReconTrackBeginningVector.Z() + fittingVectorIncrement)){
		//			if(testReconTrackBeginningVector.Mag() > radiusOfDetector) continue;
		//			for(testReconTrackEndVector.SetX(testReconTrackEndVectorStartingPosition); testReconTrackEndVector.X() <= radiusOfDetector; testReconTrackEndVector.SetX(testReconTrackEndVector.X() + fittingVectorIncrement)){
		//				for(testReconTrackEndVector.SetY(testReconTrackEndVectorStartingPosition); testReconTrackEndVector.Y() <= radiusOfDetector; testReconTrackEndVector.SetY(testReconTrackEndVector.Y() + fittingVectorIncrement)){
		//					for(testReconTrackEndVector.SetZ(testReconTrackEndVectorStartingPosition); testReconTrackEndVector.Z() <= radiusOfDetector; testReconTrackEndVector.SetZ(testReconTrackEndVector.Z() + fittingVectorIncrement)){
		//						if(testReconTrackEndVector.Mag() > radiusOfDetector || testReconTrackEndVector == testReconTrackBeginningVector) continue;

		//for(Int_t j = 0; j < plotTree->GetEntries(); j++){
		//	plotTree->GetEntry(j);
		//	testReconTrackEndVector.SetXYZ(plotVectorX, plotVectorY, plotVectorZ);
		//	if(FOM < FOMThreshold || testReconTrackEndVector == trackBeginningVector) continue;
		
		TVector3 testReconTrackDirectionVector(0., 0., 1.);		
		for(Double_t phi = 0.; phi < 2.*TMath::Pi(); phi += 0.1){
			testReconTrackDirectionVector.SetPhi(phi);
			for(Double_t theta = 0.; theta < TMath::Pi(); theta += 0.1){
				testReconTrackDirectionVector.SetTheta(theta);
								testReconTrackEndVector = trackBeginningVector + testReconTrackDirectionVector;	
								TVector3 testReconTrackUnitVector = testReconTrackEndVector - trackBeginningVector;
								testReconTrackUnitVector = testReconTrackUnitVector.Unit();
								Double_t projectionOfMaxFOM, projectionOfDeviationVector, averageOfProjections = 0.;
								Double_t deviation = 0.;
								Int_t numberOfProjections = 0;
								Double_t MaxFOM = TMath::SignalingNaN();
								for(Int_t l = 0; l < plotTree->GetEntries(); l++){
									plotTree->GetEntry(l);
									if(FOM < FOMThreshold) continue;
									TVector3 perpDeviationVector, deviationVector(plotVectorX, plotVectorY, plotVectorZ);
									deviationVector = deviationVector - trackBeginningVector;
									perpDeviationVector = testReconTrackUnitVector.Cross(deviationVector);
									deviation += FOM*perpDeviationVector.Mag2();
									projectionOfDeviationVector = testReconTrackUnitVector.Dot(deviationVector);
									averageOfProjections += projectionOfDeviationVector;
									numberOfProjections++;
									if(TMath::IsNaN(MaxFOM)){
										MaxFOM = FOM;
										projectionOfMaxFOM = projectionOfDeviationVector;
									}
									else if(MaxFOM < FOM){
										MaxFOM = FOM;
										projectionOfMaxFOM = projectionOfDeviationVector;
									}
								}
								averageOfProjections = averageOfProjections/numberOfProjections;
								if(TMath::IsNaN(deviationMin)){
									deviationMin = deviation;
									reconstructedTrackBeginningVector = trackBeginningVector;
									reconstructedTrackEndVector = testReconTrackEndVector;
								}
								else if(deviationMin > deviation){
								   	deviationMin = deviation;
									reconstructedTrackBeginningVector = trackBeginningVector;
									reconstructedTrackEndVector = testReconTrackEndVector;
								}
								if(projectionOfMaxFOM > averageOfProjections){
									TVector3 tempVector;
									tempVector = reconstructedTrackBeginningVector;
									reconstructedTrackBeginningVector = reconstructedTrackEndVector;
									reconstructedTrackEndVector = tempVector;
								}
			}
		}
		//					}
		//				}
		//			}
		//		}
		//	}
		//}

		TVector3 reconstructedTrackUnitVector = reconstructedTrackEndVector - reconstructedTrackBeginningVector;
		reconstructedTrackUnitVector = reconstructedTrackUnitVector.Unit();
//
//		//TVector3 reconstructedTrackUnitVector = trackEndVector - trackBeginningVector;
		cout
			<< "Reconstructed track direction = ("
			<< reconstructedTrackUnitVector.X()
			<< ", "
			<< reconstructedTrackUnitVector.Y()
			<< ", "
			<< reconstructedTrackUnitVector.Z()
			<< ")\n";
		TVector3 KLTrackDirection(dirnArr[0], dirnArr[1], dirnArr[2]);
		reconThetaFromRealAngle = KLTrackDirection.Angle(reconstructedTrackUnitVector);
		reconstructedTrackVectorDivergenceAngleTree->Fill();
		cout
			<< "Angle between real and reconstructed track direction vectors: "
			<< reconThetaFromRealAngle
			<< endl;
	}

	//cout << "Done reading tree." << endl;
	//dataFIn->Close();

//		TString canvasTitle;
//		ostringstream trackVectorFit;
//		trackVectorFit
//			<< "Fitted Track Vector = (" 
//			<< reconstructedTrackUnitVector.X() << ", "
//			<< reconstructedTrackUnitVector.Y() << ", "
//			<< reconstructedTrackUnitVector.Z() << "), "
//			<< "Angle difference between real and recon track = "
//			<< reconThetaFromRealAngle;
	//	cout
	//		<< "Min recon emission FOM: "
	//		<< FOMThreshold
	//		<< endl;
	//	canvasTitle = trackVectorFit.str();
	//	TCanvas* c0 = new TCanvas("c0", canvasTitle, 1200, 600);
		//TCanvas* c0 = new TCanvas("c0", "canvas 0", 600, 600);
		//c0->cd();
	//	c0->Divide(2, 1);
	//	c0->cd(1);
	//	Int_t numberOfBins = 170;
	//	Double_t detectorEdge = 850.;
		//plotTree->SetMarkerStyle(21);
		//plotTree->SetMarkerSize(0.5);
	//	plotTree->Draw("plotVectorZ:plotVectorY:plotVectorX:FOM", "FOM > FOMThreshold", "colz");
		//plotTree->Draw("plotVectorZ:plotVectorY:FOM", "", "colz");
	//	gPad->SetTheta(0.);
	//	gPad->SetPhi(0.);
	//	gPad->Update();
	//	c0->cd(2);
	//	plotTree->Draw("plotVectorZ:plotVectorY:plotVectorX:FOM", "FOM > 2000.", "colz");
	//	gPad->SetTheta(0.);
	//	gPad->SetPhi(-90.);
	//	gPad->Update();
	//	TString eventImageReconstructionFileName = "Event_image_reconstruction_with_sigma";
	//	ostringstream sigmaStream;
	//	sigmaStream << sigma;
	//	eventImageReconstructionFileName += sigmaStream.str();
	//	eventImageReconstructionFileName += ".pdf";
	//	c0->SaveAs(eventImageReconstructionFileName);
		//c0->SaveAs("positron_simulation_image.pdf");

	//	Double_t pi = TMath::Pi();
//	TCanvas* c1 = new TCanvas("c1", "canvas 1", 600, 600);
//	c1->cd();
//	reconstructedTrackVectorDivergenceAngleTree->Draw("reconThetaFromRealAngle>>thetaHistogram(20, 0., 3.14)");
//	//reconstructedTrackVectorDivergenceAngleTree->Draw("reconThetaFromRealAngle");
//	TH1D *thetaHistogram = (TH1D*)gDirectory->Get("thetaHistogram");
//	thetaHistogram->SetTitle("Deviation of Reconstructed Track Direction");	
//	thetaHistogram->SetXTitle("Deviation of Reconstructed Direction [Rad]");
//	thetaHistogram->SetYTitle("Count");
//	thetaHistogram->SetStats(1);
////	TPaveStats* histStat = (TPaveStats*)thetaHistogram->FindObject("stats");
////	histStat->SetName("reconstruction histogram");
//	c1->SaveAs("recon_track_deviation.pdf");
	cout
		<< "Program was terminated successfully!"
		<< endl;
	return 0;
}
